Transport theory is used to study the non-equilibrium physical process at the level of kinetic theory, which is crucial to study radiation transport, reactor physics, aerodynamics, plasma dynamic transport, etc. There are three major ways to descript transport theory: the micro-level description, kinetic theory description, and hydrodynamics level description. According to different physical processes, transport equations can be written in different forms, such as neutron transport equations, radiation transport equations, molecule transport equations, energetic charged particle transport equations, and plasma transport equations.

Plasma transport is widely used in astrophysics, space physics, nuclear fusion physics, etc. One of the most important aims of studying the plasma transport is to utilize the fusion energy. Plasma collective transport is a very important kind of transport in plasma physics, which is used to study the physical problems caused by collisions of two or more groups of charged particles. The basic plasma transport equations consist of continuity equation, momentum equation, and energy equation. Density control and fuel retention are two critical issues for future fusion devices with long pulse and high performance of plasma discharge.^[1] Active fuelling is a useful method to maintain the plasma density. So it is critical to improve the plasma fuelling efficiency and injection depth to meet high-performance steady-state operation requirements for the next generation of magnetic confinement fusion device. Gas puffing (GP),^[2] pellet injection (PI),^[3] and SMBI^[4] are the three major methods of plasma fuelling. Compared with GP, SMBI has a very high injection efficiency and great injection depth. Molecules can penetrate several centimeters inside the last closed (magnetic) flux surface (LCFS) during SMBI. The fuelling efficiency of SMBI is about 30%–60%. PI has a much higher fuelling efficiency, but costs much more than GP and SMBI, because of more complex techniques. SMBI can be used to reduce particle recycling^[5] and the power thresholds of L–H transition, control plasma density,^[6] mitigate the ELMs^[7,8] and study nonlocal heat transport,^[4] impurity transport,^[9] etc.

SMBI was first applied to the HL-1 tokamak device, and then installed in HL-1M,^[10] HL-2A,^[4,11] and widely in other fusion devices, such as HT-7,^[12] KSTAR,^[8,13] W7-AS,^[14] ASDEX Upgrade,^[15] JT-60U,^[6] Tore Supra,^[16] NSTX,^[17] Heliotron J,^[18] J-TEXT,^[19] and EAST.^[1] For future fusion devices especially for ITER tokamak, SMBI will also be an optional fuelling method besides PI, if its penetration depth and fueling efficiency are improved in H-mode plasma discharge. Thus, it is urgent and necessary to develop more physical models and codes of SMBI. There are many neutral transport models of fuelling such as the convection-ablation model of PI^[20–22] and the diffusion model of GP,^[23–25] but there are few models and codes of SMBI. Wang^[26] et al. developed a physical model of SMBI, which can be used to simulate the SMBI fuelling process and the GP fuelling process when the injection velocity of molecules is very low. Zhou^[27] et al. used this model to make comparisons of injection efficiency and depth between tokamak fuelling of GP and SMBI in two-dimensional (2D) simulations. Though Wang et al. have done some research about the transport of SMBI, there is still a lack of research about SMBI transport dynamics and few codes used to simulate the process of SMBI fuelling in the field. Currently, the trans_neut module of the BOUT++ code is the only large-scale parallel simulation code of SMBI in the field, while it costs lots of computer resources to run. Thousands of CPU* hours are needed to simulate a whole fuelling process with the BOUT++ code, so the efficiency of large-scale parallel simulation codes is very low. Since the radial transport dominates the SMBI transport process and determines the injection depth and efficiency, TPSMBI (TPSMBI, transport of supersonic molecule beam injection), developed in cylindrical coordinate, completely has the ability to simulate the main transport during SMBI and much higher simulation efficiency. It only costs far less than one CPU* hour to simulate a whole fuelling process. It will be an effective and necessary complimentary code for the BOUT++ code and is necessary to improve the speed of studying SMBI fuelling transport.

In this work, a new 1D code, TPSMBI, is developed completely in Fortran language, verified to be correct with some analytical solutions, and also benchmarked well with the large-scale parallel simulation code, the trans-neut module of BOUT++ code which is developed in different programing languages (such as C++) and numerical methods. TPSMBI code has many advantages: (i) it has a very high efficiency in simulation with much less computation resources; (ii) it has already captured all the key features of SMBI transport in the radial direction which is crucial to the fueling efficiency and the penetration depth; (iii) it is an effective and necessary complement to the large-scale parallel simulation code.

In this paper, we report the simulation results of molecule, atom and plasma transport during fuelling of SMBI. A specified local constant flux boundary condition is applied to model SMBI. The rest of this paper is organized as follows. The physical model is described in Section 2. The numerical method and implementation of boundary conditions are explained in Section 3. The numerical verification is shown in Section 4. The simulation results in cylindrical and slab coordinates are presented in Section 5. The summary and conclusions drawn from the present study are given in Section 6.

In tokamak plasma fuelling, hydrogen H₂ or deuterium D₂ is injected in gaseous state by GP or SMBI while it is injected in solid state by PI. During the SMBI, the molecule H₂ or D₂ will first dissociate to become H or D atoms, and then the atoms are ionized to become the major plasma. In tokamak edge plasmas, there are many kinds of collision reactions between molecules, atoms and plasmas during SMBI, which cannot be included in a simple transport model of SMBI. Only dominant reactions, such as molecule dissociation, atom ionization, and ion–atom charge exchange are included in order to reduce the difficulty in modeling the particle transport and programming TPSMBI code. This simplified model of SMBI includes four kinds of particle species (i.e., hydrogen molecules, atoms, ions, and electrons). A simple six-field fluid model, including plasma density transport equation, electron and ion heat transport equations, molecule density and momentum transport equations, and atom density transport equation, is obtained by reducing the Braginskii equations with source and sink terms. The transport equations of molecule and atom density are solved separately, since they have different kinds of sources and transport processes during SMBI. This reduced model can be used to study basic transport physics in radial direction during SMBI, such as atom density diffusion, plasma density and heat diffusion, and molecule density convection, and it captures all the key features of transport in the radial direction during SMBI. Though the model is simple, more physical reactions can be easily incorporated into the model.

Firstly, the transportation equations of plasma density N_e, ion temperature T_i, and electron temperatures T_e are studied. Plasma quasi-neutral condition (i.e., N_i = N_e) is used in the model. The equations are as follows:

In the molecule dissociation process and atom ionization process, electrons provide energy (i.e., W_diss) to divide molecules into atoms and ionize atoms. In the charge change process, the energy exchange takes place between ions and atoms.

Secondly, we consider the atom density (N_a) transport equation, on the basis of atom diffusion caused by the strong charge exchange collision rate in the model. The atom and ion temperatures are typically assumed to be equal, i.e. T_i = T_a. The atom density transport equation is

Thirdly, molecule transport equations consist of molecule density N_m transport equation and radial velocity V_rm transport equation. The equations are as follows:

Definitions of the sources and rates due to particle collision reactions are as follows:

The newly developed code TPSMBI is used to simulate the SMBI fuelling process and other transport processes.

3.1. Numerical methods in TPSMBI

The convection–diffusion equations solved by TPSMBI have the following simplified form

where u is the physical field, V is the convection velocity, D is the diffusion coefficient, and S is the source or sink of the physical field. The time difference scheme is the θ scheme in TPSMBI. The second-order central difference scheme for the diffusion term and first-order upwind scheme for the convection term are adopted. Besides, other numerical methods for convection term, such as the second-order upwind method, third-order upwind method, are optional to improve numerical precision if necessary in TPSMBI code. The upwind scheme is stable if the following Courant–Friedrichs–Lewy (CFL) condition is satisfied:

where DT is the time interval, and DX is the space interval. Modified wavenumber analysis shows that the first-order upwind scheme introduces severe numerical diffusion in the solution where large gradients exist. The second-order upwind scheme is less diffusive than the first-order accurate scheme, while it introduces dispersive errors in the region where the gradient is high. The third-order upwind scheme is less diffusive and dispersive than the second-order scheme. Thus, the third-order upwind scheme is used in simulation. In TPSMBI, spatial difference terms are divided into two parts, and their weights are θ and (1 − θ), respectively. The method in TPSMBI code is unconditionally stable if θ ≥ 1/2. For θ = 1/2, a Taylor series analysis of the scheme (Crank–Nicolson scheme) shows that it is second-order precision in space and time.

There are finite difference steps in space and time. Thus denote x_j = x₀ + j * DX, j = 0,1,…,J and t_i = t₀ + i* DT, i = 0,1,…,I. In Eq. (14), applying the first order upwind method to the convection term and second-order central difference method to the diffusion term, we have

Equation (16) can easily be rewritten as follows:

where and .

It is convenient to rewrite Eq. (17) in matrix form as

where A is a tri-diagonal matrix as follows:

Equation (18) is solved with Gauss elimination with a partial pivoting method, which has much higher precision.

3.2. Numerical boundary conditions

In simulations, we choose N_m0 = N₀ = 1 × 10¹⁹ m⁻³ and V_rm0 = − 1000 m/s according to the experiment parameters of SMBI in HL-2A as the edge boundaries of the molecule density and the radial molecule velocity. Neumann boundary conditions are applied to the core boundary of molecule density and radial velocity. At the core and edge boundaries, atom density is given in Neumann boundary conditions. The boundary conditions of plasma density, electron and ion temperatures are as follows:

With these conditions, stable profiles of plasma density, electron and ion temperatures can be obtained, which are consistent with the results of tokamak experiments. Because the slopes of the profiles of electron density, electron and ion temperatures at steady state are the same as and , respectively, in 1D slab coordinate. Thus, it is convenient and efficient for TPSMBI to initialize the profiles of electron density, electron and ion temperatures, which are the same as their stable profiles determined by the boundaries.

Transport dynamics during SMBI is studied after the system has reached its steady state, which is set as t = 0.00 ms in simulation. SMBI is used during the whole simulation period with a constant molecule density and an inward constant radial velocity V_rm0 = − 1000 m/s at the last closed flux surface.

Firstly, transport dynamics of molecule density and atom density are studied in the radial direction during SMBI. Inward propagations of molecules and atoms are observed during SMBI in both cylindrical coordinate and slab coordinate as shown in Figs. 5 and 6. There are neither molecules nor atoms at t = 0.00 ms. The time intervals between two solid (dash) curves are the same in Fig. 5, and the distance between two curves decreases, which means the propagation speed of the molecule front becomes smaller as molecules propagate into the plasma more deeply due to the increased molecule dissociation rate. The molecule density at the leading edge of the front in the cylindrical coordinate at t = 0.05 ms is smaller than at the source boundary due to the source of molecule dissociation. While the molecule densities at the leading edge of the molecule front at t = 0.15 ms and t = 0.20 ms are larger than at the boundary and the difference increases with molecule beam propagating into the plasma due to the effect of the term (1/r)V_rmN_m in Eq. (5). This term only appears in the molecule density equation in the cylindrical coordinate but not in the slab coordinate. It acts like a source in the molecule density equation in the cylindrical coordinate and leads to the different simulation results between slab and cylindrical coordinates as shown in Fig. 5. The penetration depth of molecule in cylindrical coordinate is deeper than in slab coordinate due to the curvature effect or the term (1/r)V_rmN_m. The difference becomes larger once it propagates more deeply with time increasing due to the increase of the curvature effect. The physical picture of the different simulation results of molecule density between slab coordinate and cylindrical coordinate can be like that along the molecule propagation path, the radius r and the ring area between r and r + Δr decrease. In order to keep the number of molecules in the reducing ring area unchanged along the propagation path, the molecule density increases. The increased molecule density makes molecule propagate more deeply in the cylindrical coordinate than in the slab coordinate. The difference in molecule density between in the slab coordinate and in the cylindrical coordinate increases as radius r decreases, and decreases as the radius r increases. If the minor radius a increase from 0.4 m to 2.05 m which means that the radius and the ring area increase, the difference between slab and cylindrical coordinates will become much smaller, and even can be neglected. Once the minor radius goes to infinity, the curvature effect in cylindrical coordinate can be neglected in the edge domain simulation (i.e., , which means that it can be treated in the way the same as that in slab coordinate.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Time evolution of radial profiles of molecule density in slab geometry (dash curves) and cylindrical geometry (solid curves).

The molecule transport is dominated by the convection due to the directed injection velocity of molecule, while the atom transport is dominated by the atom diffusion. The increase of atom density at edge and propagating depth slightly deeper than that of molecule can be due to the diffusion of atoms inside the last closed flux face. The reason why the atom density is larger in the cylindrical coordinate than in the slab coordinate can be mainly due to the curvature effects and the larger molecule density in the cylindrical coordinate than in the slab coordinate as shown in Fig. 6. So is the electron density as shown in Fig. 7.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Time evolution of radial profiles of atom density in cylindrical coordinate (solid curves) and in slab coordinate (dash curves).

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Time evolution of radial profiles of electron (plasma) density in 1D cylindrical coordinate (solid curves) and in slab coordinate (dash curves).

In the inward propagation period of molecules and atoms, plasma density increases locally, while electron and ion temperatures decrease locally along the molecule propagation path as shown in Figs. 7–9. Plasma density increases dramatically near the beam front as shown in Fig. 7, due to ionization of atoms and the much smaller diffusion coefficient than atoms’. The decrease of electron temperature can be mainly due to electron energy loss in the processes of dissociation and ionization, while ion temperature decreases due to the increase of ion density and energy interchange between electrons and ions. The electron density is larger in the cylindrical coordinate than in the slab coordinate. The electron and ion temperatures are slightly lower in the cylindrical coordinate than in the slab coordinate.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Time evolution of radial profiles of electron temperature in 1D cylindrical coordinate (solid curves) and in slab coordinate (dash curves).

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. Time evolution of radial profiles of ion temperature in 1D cylindrical coordinate (solid curves) and in slab coordinate (dash curves).

A plasma fuelling model of SMBI including molecules, atoms and plasmas is developed in 1D cylindrical coordinate. The dominant collisional interactions between molecules, atoms and plasmas, such as molecule dissociation, atom ionization and atom–ion charge exchange are included in the model. A new code, TPSMBI code, is developed completely, verified to be correct with some analytical solutions and also benchmarked well with a larger-scale parallel trans_neut module of BOUT++ code. Radial transport dynamics of SMBI is simulated with TPSMBI code in both slab and cylindrical coordinates. SMBI is modeled by giving a constant molecule flux boundary condition at the last closed flux surface. Flux-driven boundary conditions are applied to plasma density, electron and ion temperatures at the core of the simulation region to maintain the profiles. Principal results are summarized as follows.

With a constant injection flux of molecule density at the last closed flux surface, molecules and atoms can propagate inward about 7 cm inside the last closed flux surface in cylindrical coordinate. The propagation speed of the front of molecules becomes smaller with increasing dissociation rate of molecule along the molecule propagation path. The plasma density keeps increasing locally and plasma temperature keeps decreasing locally due to particle fuelling, heat sink, and diffusion during SMBI.

The different simulation results obtained in cylindrical and slab coordinates are mainly due to the different molecule densities caused by the curvature effect and the reducing ring area between r and r + Δr along the molecule propagation path in cylindrical coordinate.